# Alexander F. Gazmararian
# afg2@princeton.edu

# Load packages
library(tidyverse)
library(scales)
library(modelsummary)
library(dfadjust)
library(emmeans)
library(margins)

# Load functions
source("code/fun/book_theme.r")
source("code/fun/savefig.r")
source("code/fun/fix_txt.r")

# Load data
g <- readRDS("data/NatRegQual_Winter21.rds")

# Subset to the credibility sample (the regulation experiment is not on the Roosevelt sample)
g <- subset(g, sample == "credibility")

# Specify baseline model
source("code/fun/model_spec.r")

# Create coef names map
source("code/fun/coefnames4tables.r")
coefnames <- c(
  "support_assist"="Bipartisan",
  "reverse_assist"="Not Reversible",
  "support_assist:reverse_assist"="Bipartisan x Not Reversible",
  "support_assist:PartySummaryNeither"="Bipartisan x Independent",
  "support_assist:PartySummaryRepublican"="Bipartisan x Republican",
  "reverse_assist:PartySummaryNeither"="Not Reversible x Independent",
  "reverse_assist:PartySummaryRepublican"="Not Reversible x Republican",
  "support_assist:reverse_assist:PartySummaryNeither"="Bipartisan x Not Reversible x Independent",
  "support_assist:reverse_assist:PartySummaryRepublican"="Bipartisan x Not Reversible x Republican",
  coefnames
)

# Create Figure 5.6 ----
m1 <- lm(reg_support_i ~ support_assist, g)
m1.mean <- emmeans(m1, ~ support_assist)
m1.mean <- data.frame(m1.mean)
m1.mean$treatment_value <- c("Low", "High")
m2 <- lm(reg_support_i ~ reverse_assist, g)
m2.mean <- emmeans(m2, ~ reverse_assist)
m2.mean <- data.frame(m2.mean)
m2.mean$treatment_value <- c("Reversible", "Not Reversible")
names(m2.mean)[1] <- "treatment"
names(m1.mean)[1] <- "treatment"
bipart <- rbind(m1.mean, m2.mean)
bipart$lb90 <- with(bipart, SE * qt(.05, df) + emmean)
bipart$ub90 <- with(bipart, SE * qt(.95, df) + emmean)
bipart$model <- c("Bipartisanship", "Bipartisanship", "Institutional Constraints", "Institutional Constraints")
bipart$treatment_value <- factor(bipart$treatment_value, ordered = TRUE, levels = c("Low", "High", "Reversible", "Not Reversible"))
plot.desc <-
  bipart %>%
  ggplot(aes(x=treatment_value,y=emmean,color=factor(treatment),label=scales::percent(emmean, accuracy=1))) +
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width = 0) +
  geom_errorbar(aes(ymin=lb90, ymax=ub90), width = 0, size = 1.75) +
  geom_point(size = 4) +
  facet_wrap(~model, scales = "free_x") +
  geom_text(hjust=-.75, size = 3, color = "black") +
  scale_color_grey() +
  book_theme +
  scale_y_continuous(label = scales::percent, limits = c(.55,.82), expand=c(0,0)) +
  labs(x="Treatment Condition", y = "Support") +
  theme(legend.position="")
plot.desc
savefig(plot.desc, "5.6_figure_regsup", height = 2, filepath = "figures/")

# Create table for online appendix ----

# Analyze distribution of treatment
prop.table(table(g$reg_treat_cat))

# Analyze distribution of response by the treatment
prop.table(table(g$reg_support, g$reg_treat_cat), 2)

# Examine results by partisanship 
with(subset(g, PartySummary=="Republican"), prop.table(table(reg_support, reg_treat_cat), 2))
with(subset(g, PartySummary=="Democrat"), prop.table(table(reg_support, reg_treat_cat), 2))
with(subset(g, PartySummary=="Neither"), prop.table(table(reg_support, reg_treat_cat), 2))

# Estimate models
m <- list()
# Binary outcome
m[[1]] <- lm(reg_support_i ~ support_assist * reverse_assist, g)
m[[2]] <- lm(update(f.base, reg_support_i ~ support_assist * reverse_assist + .), g)
m[[3]] <- lm(update(f.base, reg_support_i ~ support_assist * reverse_assist + . + gw_index_cred + pc1), g)
m[[4]] <- lm(update(f.base, reg_support_i ~ support_assist * reverse_assist * PartySummary + . + gw_index_cred + pc1), g)
# Scale outcome
m[[5]] <- lm(as.numeric(reg_support) ~ support_assist * reverse_assist, g)
m[[6]] <- lm(update(f.base, as.numeric(reg_support) ~ support_assist * reverse_assist + .), g)
m[[7]] <- lm(update(f.base, as.numeric(reg_support) ~ support_assist * reverse_assist + . + gw_index_cred + pc1), g)
m[[8]] <- lm(update(f.base, as.numeric(reg_support) ~ support_assist * reverse_assist * PartySummary + . + gw_index_cred + pc1), g)

# Create table
file <- "tables/ch5/ols_reg.txt"
modelsummary(
  m,
  stars = c("*"=.1,"**"=.05,"***"=.01),
  vcov = "HC2",
  gof_map = c("nobs","adj.r.squared"),
  coef_map = coefnames,
  output = "latex",
  escape = FALSE
) %>%
  cat(., file = file)
fix_txt(file)

# Create online appendix table for interpretation check ----

# Inspect outcome variable
with(g, prop.table(table(reg_reverse, reg_treat_cat), 2))

# Estimate models
m <- list()
# Binary outcome
m[[1]] <- lm(reg_reverse_i ~ support_assist * reverse_assist, g)
m[[2]] <- lm(update(f.base, reg_reverse_i ~ support_assist * reverse_assist + .), g)
m[[3]] <- lm(update(f.base, reg_reverse_i ~ support_assist * reverse_assist + . + gw_index_cred + pc1), g)
m[[4]] <- lm(update(f.base, reg_reverse_i ~ support_assist * reverse_assist * PartySummary + . + gw_index_cred + pc1), g)
# Scale outcome
m[[5]] <- lm(as.numeric(reg_reverse) ~ support_assist * reverse_assist, g)
m[[6]] <- lm(update(f.base, as.numeric(reg_reverse) ~ support_assist * reverse_assist + .), g)
m[[7]] <- lm(update(f.base, as.numeric(reg_reverse) ~ support_assist * reverse_assist + . + gw_index_cred + pc1), g)
m[[8]] <- lm(update(f.base, as.numeric(reg_reverse) ~ support_assist * reverse_assist * PartySummary + . + gw_index_cred + pc1), g)

# Create table
file <- "tables/ch5/ols_reg_interp.txt"
modelsummary(
  m,
  stars = c("*"=.1,"**"=.05,"***"=.01),
  vcov = "HC2",
  gof_map = c("nobs","adj.r.squared"),
  coef_map = coefnames,
  output = "latex",
  escape = FALSE
) %>%
  cat(., file = file)
fix_txt(file)
